#!/bin/csh -f
#goal: read CFAD.bin    
#      then draw vertical profile CFAD picture
#  by: wangh 2015/3/18
#-----------------------------------------
#set debug
set cmd = $0
set cmd = $cmd:t
#-----------------------------------------
if ( $#argv == 1 ) goto usage

switch ( $#argv )
   case 0:
   breaksw
  default:
usage:
  echo " "
  echo "Usage : $cmd no input parameter"
  echo " "
  exit(2)
  breaksw
endsw
#-----------------------------------------
set ExpNam = "FR"
set HomeDir = "/public/home/wangh29/data2019-2021/"
set wrfNC  = "/public/home/wangh29/data/model/THOM/wrfout_d03_2017-05-06_17_00_00.nc"
set outdir  = "./svg/fig15.SketchMap${ExpNam}"
set SecLon = 112.5
set latS   = 20.
set latE   = 24.5

set RadarS = "TRUE"
set TempL  = "TRUE"
set CoastL = "TRUE"
#-----------------------------------------
#-----------------------------------------
#Parameter set 
#-----------------------------------------
set type = "x11" #pdf,x11,ps,ncgm
set FontS = 0.025
#-----------------------------------------
#Create plot.ncl
#-----------------------------------------
  cat >obs_radar.ncl << EOF
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/panel_two_sets.ncl"
;  These files still have to be loaded manually
;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
  begin
; We generate plots, but what kind do we prefer?
;-------------------------------------------
   wks_type ="${type}"
   if("${type}".eq."png")then
    wks_type@wkWidth = 2500
    wks_type@wkHeight = 2500
   else
    wks_type@wkPaperHeightF  =11 
    wks_type@wkPaperWidthF = 11 
   end if 
   wks  = gsn_open_wks(wks_type,"${outdir}")
   gsn_define_colormap(wks,"NCV_jet")
;-------------------------------------------
   a = addfile("${wrfNC}","r")
   zs=  wrf_user_getvar(a,"HGT",0)
   lon  = wrf_user_getvar(a,"XLONG",0) 
   lat  = wrf_user_getvar(a,"XLAT",0)
   lon1d = dim_avg_n(lon,0)
   lat1d = dim_avg_n(lat,1)
   SecID = toint((${SecLon}-lon1d(0))/0.01)
   print("section at "+lon1d(SecID))
   zs@_FillValue = -999.
   zs= smth9(zs,0.50,  0.25, True)
   zs1d = dim_avg_n(zs(:,SecID-100:SecID+100),1)/1000.
   ns = dimsizes(zs1d)
   zs2d = new((/2,ns/),"float")
   zs2d(0,:) = 0.
   zs2d(1,:) = zs1d
;----------------------------------------
   if("${ExpNam}".eq."FR")then
    LeftStr = " a) ${ExpNam}"
   elseif("${ExpNam}".eq."WR")then
    LeftStr = " b) ${ExpNam}"
   end if 

  res = True
;  res@vpWidthF = .6
;  res@vpHeightF= .5
  res@gsnMaximize = True
  res@gsnDraw     = False
  res@gsnFrame    = False
  res@gsnLeftString = LeftStr;
  res@gsnLeftStringOrthogonalPosF = -0.1
  res@gsnCenterString  = "  "
  res@gsnRightString  = "  "
  res@gsnRightStringOrthogonalPosF = -0.1

  res@gsnLeftStringFontHeightF   = ${FontS} 
  res@gsnRightStringFontHeightF  = ${FontS} 
  res@gsnCenterStringFontHeightF = ${FontS}
  res@tiYAxisOn             = True;False 
  res@tiYAxisString         = "Height (km)"
  res@tiXAxisOn             = True 
  res@tiXAxisString         = "Lattitude" 
  res@tiXAxisFontHeightF    = ${FontS}
 
  res@trYMinF = 0.
  res@trYMaxF = 13.
  res@trXMinF = 19.
  res@trXMaxF = 25.9

 ;tickmark res
  res@tmXTOn                 = True;False
  res@tmXTLabelsOn = False
  res@tmXBLabelsOn = True
  res@tmXBLabelFontHeightF  = ${FontS};0.015
  res@tmYLLabelFontHeightF  = ${FontS};0.015
  res@tmYLLabelStride = 2
  res@tmXBLabelStride = 2

  res@gsnXYFillColors = "brown";"dim gray"

   res@xyLineThicknesses   = (/3.0/)              ; make line thicker
   res@xyLineColors        = (/"transparent","brown"/);           ; change line color
   plot  = gsn_csm_xy(wks,lat1d,zs2d,res) 

;-----------------------------------------------------
; precipitation area
  poly_res = True
  poly_res@gsFillColor = "dim gray"    ; set fill color to white
  if("${ExpNam}".eq."FR")then
	  gsn_polygon_ndc(wks,(/0.57,0.77,0.77,0.57,0.57/),(/0.15,0.15,0.175,0.175,0.15/),poly_res)
  elseif("${ExpNam}".eq."WR")then
  gsn_polygon_ndc(wks,(/0.43,0.59,0.59,0.43,0.43/),(/0.15,0.15,0.175,0.175,0.15/),poly_res)
  end if 
;--------------------------------
;text and marker resourse
  txres = True
  txres@txFontHeightF = ${FontS}*0.9
  txres@txFontColor  = "black"
  pmres = True
  pmres@gsMarkerSizeF = 0.03
  pmres@gsMarkerOpacityF = 1.0
  pmres@gsMarkerColor = "black"
  pmres@gsMarkerIndex = 12
  pmres@gsMarkerThicknessF = 2 
;--------------------------------
;radar line 
 if("${RadarS}".eq."TRUE")then 
  dpr_lat = (/23.0039,21.845/)
  dpr_hgt = (/0.2,0.2/)
  dpr_name=(/"GZ","YJ"/)
  MarkGZ= gsn_add_polymarker(wks,plot,dpr_lat,dpr_hgt,pmres)
  textGZ= gsn_add_text(wks,plot,dpr_name,dpr_lat+0.34,0.2,txres)
 end if 
;-----------------------------------
;coastline mark
 if("${CoastL}".eq."TRUE")then
  CoastLat = (/21.445/)
  pmres@gsMarkerIndex = 8
  pmres@gsMarkerSizeF = 0.02
  pmres@gsMarkerThicknessF = 4
  pmres@gsMarkerColor = "blue"
  MarkCoast= gsn_add_polymarker(wks,plot,CoastLat,0.2,pmres)
  txres@txAngleF = 45.
  txres@txFontColor= pmres@gsMarkerColor
  textGZ= gsn_add_text(wks,plot,"Coastline",CoastLat-0.34,0.2,txres)
 end if 
;------------------------------------
;line resourse
  AvgFrlv = 5.11
  Avg20Hgt= 8.85
  resp                  = True              ; polyline mods desired
  resp@gsLineColor      = "dim gray"           ; color of lines
  resp@gsLineThicknessF = 4.0               ; thickness of lines
  resp@gsLineDashPattern= 0
  if("${TempL}".eq."TRUE")then
  xpts = (/-100., 200./)
  ypts = (/AvgFrlv,AvgFrlv/)
  ypts2 = (/Avg20Hgt,Avg20Hgt/)
  dumL  = new((/2/),"graphic")
  resp@gsLineDashPattern= 0
  dumL(0) = gsn_add_polyline(wks,plot,xpts,ypts,resp)
  resp@gsLineDashPattern= 2
  dumL(1) = gsn_add_polyline(wks,plot,xpts,ypts2,resp)
  txres@txAngleF = 0.
  txres@txFontColor= resp@gsLineColor
  text0= gsn_add_text(wks,plot,"0~S~o~N~C",19.5,AvgFrlv+0.6,txres)
  text2= gsn_add_text(wks,plot,"-20~S~o~N~C",19.5,Avg20Hgt+0.6,txres)
  end if 


   draw(plot)
   frame(wks)


end
EOF

ncl obs_radar.ncl 
rm  obs_radar.ncl

